#include "cutrace.h"
#include <math.h>
#include <cstdlib>
#include <stdio.h>

extern double Bfr[Nx*Ny];
extern double Bfz[Nx*Ny];
extern double Bfphi[Nx*Ny];
extern double Apsi[Nx*Ny];

void initdata(vec* vap, vec* vav, int N)
{
    double energymax = 0.5;  
    int Nr = 1 << 3;
    int Nz = 1 << 3;
    int Nvpar = 1 << 3;
    int Nvper1 = 1 << 3;
    int Nvper2 = 1 << 3;
    for (int i = 0; i < N; i++)
    {
        double ri = 0.4*(i&7)/(double)Nr;
        double zi = 0.4*((i&(7<<3))>>3)/(double)Nz;
        double vpar = energymax*((i&(7<<6))>>6)/(double)Nvpar;
        double vper1 = energymax*((i&(7<<9))>>9)/(double)Nvper1;
        double vper2 = energymax*((i&(7<<9))>>9)/(double)Nvper2;
        vap[i] = {ri+3.8, 0.0, -0.2+zi};
        vec vB = vectorfield(vap[i]);
        vec vb = 1.0/pow(dot(vB,vB),0.5)*vB;
        vec a = {0.0,0.0,1.0};
        vec e1 = 1.0/pow(dot(cross(vb,a),cross(vb,a)),0.5)*cross(vb, a);
        vec e2 = cross(vb,e1);
        vav[i] = vpar*vb + vper1*e1 + vper2*e2;   
    }
    return;
}

void initdata_passing(vec* vap, vec* vav, int N)
{
    double dx = 0.8/(double)N;
    double R0 = (Rmax+Rmin)/2.0;
    for (int i = 0; i < N; i++)
    {
        vap[i] = {R0+dx*i, 0.0, 0.0};
        vav[i] = {0.0,0.5,0.3};
        if(checkboundary(vap[i]) == false)
        {
            printf("particle outside boundary\n");
            continue;
        }
    }
    return;
}

double func(double E, double pitch, double psi)
{
    double f = 1.0;
    double v = sqrt(E);
    double vcut = 1.2;
    double vcrit = 0.5;
    double pitch0 = 0.8;

    f = f*1.0/(pow(v,3)+pow(vcrit,3));
    f = f*exp(-pow((pitch-pitch0),2));
    f = f*exp(-psi/0.25);

    return f;
}  

void initdata_distribfunc(vec* vap, vec* vav, int N)
{
    srand(23);
    for(int i = 0; i < N;)
    {
        double R = Rmin + (Rmax-Rmin)*rand()/512;
        double Z = Zmin + (Zmax-Zmin)*rand()/512;
        double v = rand()/512;
        double psi = interp(R, Z, Apsi);
        double prob = func(v*v, 0.8, psi);
        if (prob > 0.8)
        {
            vap[i] = {R, 0.0, Z};
            vav[i] = {v, v, v};
            i = i+1;
        }
    }
}

